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The first-principles description of liquid water using ab initio molecular dynamics (AIMD) 
based on Density Functional theory (DFT) has recently been found to require long equili- 
bration times, giving too low diffusivities and a clear over-structuring of the liquid. In the 
light of these findings we compare here the room-temperature description offered by two dif- 
ferent exchange correlation functionals: BLYP, the most popular for liquid water so far, and 
RPBE, a revision of the widely used FEE. We find for RPBE a less structured liquid with 
radial distribution functions closer to the experimental ones than the ones of BLYP. The 
diffusivity obtained with RPBE for heavy water is still 20% lower than the corresponding 
experimental value, but it represents a substantial improvement on the BLYP value, one 
order of magnitude lower than experiment. These characteristics and the hydrogen-bond 
(HB) network imperfection point to an effective temperature ^^3% lower than the actual 
simulation temperature for the RPBE liquid, as compared with BLYP's ^17% deviation. 
The too long 0-0 average nearest-neighbor distance observed points to an excessively weak 
HB, possibly compensating more fundamental errors in the DFT description. 

I. INTRODUCTION 

The interaction of organic pollutant molecules with the minerals and other components of the 
soil represents a paradigmatic case within the set of atomic-scale problems of relevance for the 
environment. It displays the level of complexity that demands a fundamental revision of the way 
we use computers for the simulation of condensed matter. On one hand the number of adsorption 
studies to address is enormous (combinations of many different pollutants, minerals, surfaces thereof 
etc.), on the other, each study displays a large degree of complexity, not least because of the wet 
character of these systems. 

The effect of water in this and other wet systems is crucial in different ways, from the dielectric 
screening of electrostatic interactions to the active chemical role played in some important wet 
processes. These different aspects of water can be described at different levels of theory, including 
the static calculation of the system of interest immersed in an effective dielectric continuum, or the 
dynamic simulation of the liquid based on empirical potentials of different kinds. Even if simplified 
descriptions of water can be appropriate for some problems, other problems require the accuracy 
and transferability provided by first-principles simulation methods. This is not only our case, but 
it is true for many wet systems throughout geochemistry, biochemistry and wet chemistry. These 
considerations and the high variety of situations that we encounter in our pollutant problem (see 
Ref. Q 

for the combinatorial aspects of the project and the workflow issues it entails) make it 
extremely important that we can use different simulation methods that involve different compu- 
tational demands in a seamless fashion and within a single computational environment We 
need not only to interconnect the different simulations^ but that each one of them performs at the 
required level of accuracy. 

On the first-principles side of the wet-systems description, ab initio molecular dynamics (AIMD) 
based on density-functional theory (DFT) is the technique that, on one hand, satisfies the needs 
mentioned above, and, on the other, is efficient enough to allow for the system sizes and time 
scales needed in the simulations with today's computers. Indeed, the last ten years have witnessed 
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the success of many DFT-based AIMD simulations of liquid water, including work on structural, 
dynamical, chemical and electronic properties i^i^i^i^'^i'^^i^'^i^^ 

Recent reportsfi^^'^^^ however, have questioned some of the results of earlier studies, showing 
that if the simulations are allowed to run for longer times, the diffusivity drops by one order of 
magnitude and the liquid becomes over-structured. The discrepancy with experiment a^7i^?i^? is still 
unexplained. Prime suspects are the fundamental limitations of present-day AIMD simulations of 
liquid water: the inability of gradient-corrected (GGA) density functionals to describe dispersion 
interactions, and the neglect of quantum fluctuations in the classical description of the nuclear 
dynamics. 

The need for long equilibration times and the poor structural and dynamical results have been 
further assessed in a recent contribution,^'^ where a relaxation process of a time scale larger than 
20 ps was found at room temperature. It was observed to be related to rearrangements of the 
hydrogen-bond (HB) network in the liquid, expressed by the concentration of under-coordinated 
molecules, including a substantial presence (and relevance) of bi-coordinated ones. Both the radial 
distribution functions (RDFs) and the diffusivity were found to correlate strongly with the HB 
network imperfection, in close analogy to what Giovambattista et found in supercooled water. 
Furthermore, while the network was slowly evolving, the diffusivity was observed to equilibrate 
much faster to the instantaneous network state than to the final temperature. The mentioned 
paper also shows that the structure, diffusivity and network characteristics obtained, would very 
nicely fit experiment if they corresponded to an effective temperature of ~240 K, instead of the 
300 K actually used in the AIMD simulations. 

All those results were obtained for gradient-corrected (GGA) density functionals, specifically 
for two proposals, BLYP^^'"^^ and PBE,^ Both very popular and quite successful, the former 
comes from the quantum-chemistry community and the latter from the condensed-matter physics 
community. Even if they were constructed with different criteria, from a fundamental point of view 
it is good news that both functionals give very similar results*^ From a more practical view point, 
however, it would be highly desirable to have a satisfactory and efficient DFT description of liquid 
water in order to be able to address complex wet systems, especially since the more technical aspects 
for efficient DFT calculations (including linear scaling) have already been successfully tested for 
this system^i 

Alternative GGAs are being explored at the moment, and very promising results have been 
found by VandeVondele et al^, who tried several functionals finding much more satisfactory re- 
sults. Since their functionals originate from the chemistry community, where special emphasis is 
put on accuracy for light elements, and since our interest in wet mineral surfaces involves heavier 
elements, we think it is worthwhile to complement their work with a similar study with functionals 
proposed more generally for the whole periodic table. In the following we present results for the 
structure, diffusivity and equilibration times for RPBE.^- This GGA functional is a modification 
of PBE, well within its philosophy and fulfilling all its fundamental criteria. It is also very ac- 
cessible from the points of view of coding (only a few lines differ from a PBE implementation) 
and computational demands, remaining at the GGA level. It has already been tested for varied 
systems, including systems with heavy elementsi^i^S Functionals beyond GGA (meta-GGAs) j^?'?? 
hybrid functionalsj^L^Si^ and functionals including Van der Waals interactionaSiiS^ should also be 
checked, but their efficiency is lower, and the scope of wet systems that can be addressed at present 
is much narrower than for GGAs. None of these functionals will be addressed here. 

Asthagiri et alM- tested an alternative GGA, very similar to PBE as well, called revPBEi^ The 
spirit of the PBE modification of this functional is similar to RPBE (except for the fact that the 
latter respects the Lieb-Oxford bound locally, which ensures a global fulfilment of the condition 
for any electronic density) representing an alternative worth considering. Good revPBE results for 
liquid water were obtained by those authorsi^ before the long time-scale relaxation problem was 
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TABLE I: AIMD simulations performed in this work, all of them with periodic boundary conditions in a 
cubic cell of size a=9.865 A. DF stands for the particular density fimctional used, T for the final equilibrated 
temperature, Tsim for the AIMD simulation time after AIMD equilibration, Tgg for the AIMD equilibration 
time, "Model" for the model used for preparation, Tpre for the temperature at which the preparation model 
had been equilibrated and Ti for the AIMD initial temperature (after the Teg anneal). Temperatures in K 
and times in ps. 
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30 
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325 


325 



reported. Under the new circumstances it seems that the equilibration and simulation times used 
therein could have been too short. A revision of those tests would also be useful, but is beyond 
the scope of the present work. 

II. METHOD 

The simulations were performed using the Kohn-Sham approach'^^ to DFT— in the generalized- 
gradient approximation (GGA). The BLY P^^i^? and RPBE^^^ exchange-correlation functionals 
were used and are compared in the following. Core electrons were replaced by norm-conserving 
pseudopotentials"^^ in their fully non-local representation.^*^ Numerical atomic orbitals (NAO) of 
finite support were used as basis set, and the calculation of the self-consistent Hamiltonian and 
overlap matrices was done using the linear-scaling Siesta methodi^i^ Integrals beyond two-body 
were performed in a discretized real-space grid, its fineness determined by an energy cutoff of 150 
Ry. A double-^ polarized (DZP) NAO basis set was used, which had been obtained following the 
method proposed in Refs. 14^13 for a confining pressure^ of 0.2 GPa. The validation of the method, 
pseudopotentials and basis set can be found in Ref. [i^, including the approximations required for 
linear scaling. 

We performed AIMD simulations for 32 molecules of heavy water for the BLYP and RPBE den- 
sity functionals (a detailed comparison between BLYP and PBE for this system can be found in 
Ref.Ei). Details are given in Tabled AIMD equilibration is achieved by means of temperature an- 
nealing (velocity re-scaling)^ while the actual simulations are performed by Verlet's integration!^ 
In all the simulations the time step used was 0.5 fs. The observed total-energy drifts corresponded 
to drifts in the system temperature between 0.26 K/ps and 0.36 K/ps. The simulations were 
performed at constant volume, i.e., for fixed cell size and shape, under periodic boundary con- 
ditions. As in Ref. H the time scale of the HB relaxation process is monitored by following its 
non-equilibrium behaviour. The "moving time window" used in Ref. [i^ is used to obtain the time 
dependence of the variables of interest. 

Empirical simulations were performed using different force fields (TIP5P^ and SPC/E^^) in 
order to prepare reasonably equilibrated starting points for AIMD. These simulations were per- 
formed with the GROMACS MD package^*^'*^^ under constant volume and temperature conditions 
using a Berendsen-type thermostat*^ The empirical simulations were equilibrated during 200 ps. 
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FIG. 1: Comparison of the H-H, 0-H, and 0-0 radial distribution functions as obtained in this work for 
32 molecules with RPBE (solid line), with BLYP (dotted line) and by experimen t i ^ i " i mm (dashed line), at a 
temperature of 300 K. 

III. RESULTS AND DISCUSSION 

Fig. ^ compares the RDFs obtained with RPBE and BLYP to experiment. Looking mainly at 
the inter- molecular features, the RDFs for RPBE are considerably less structured than the ones for 
BLYP, even if still over-structured with respect to experiment. The comparison in Ref. 03 between 
PBE and BLYP gave much closer results. The too large distance for the first peak position of the 
0-0 RDF is worth noting for RPBE, however. It indicates too weak HBs that could account for 
the more fluid character of the RPBE liquid as compared with BLYP or PBE. The fact that this 
0-0 distance is longer than the experimental value (possibly weaker HB) but the liquid remains 
still slightly overstructured otherwise, suggests the presence of a more fundamental error (possibly 
the same as for BLYP and PBE) that is partly compensated by this weaker HB. 

The values for the corresponding diffusivities are shown in Table ITll The table also shows the 
effective temperatures^^ for the different simulations, i.e., the temperature at which the diffusivity 
of the real liquid equals the AIMD value. The table shows the effective temperature with respect 
to both normal and heavy water. In Ref. ii20 there was no need to refer to the difference between 
heavy and light water, since the differences in diffusivity between the two (~25%) are negligible 
when compared with the deviation of AIMD (BLYP) with respect to experimenti^ In this study, 
however, the RPBE description is much closer to experiment, and the distinction becomes relevant. 

As shown in Ref.Eol. several results for the room-temperature simulations fairly reproduce results 
of experiments or empirical models at the effective temperature. Besides the diffusivity, features 
of the RDFs like the heights of the first minimum and the second maximum of the 0-0 RDF 
would conform to this temperature scaling. The HB network also follows this scaling. In Fig. IJJthe 
distribution of molecular coordinations is compared for RPBE and TIP5P at room temperature and 
BLYP at 345 K. The RPBE histogram shows a very slightly more over-coordinated liquid than the 
TIP5P one, suggesting a slightly lower effective temperature than 300 K. A remarkable agreement 
appears when comparing the diffusivities of both simulations, BLYP at 345 K and RPBE at 300 K. 
As reported in Ref. the diffusivity of the former is ~ 1.5 x 10~^ cm^/s, identical to the one 
found for the latter in this work. 

Fig.Elshows the evolution of the diffusivity and relates it to the evolution of the liquid structure. 
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TABLE II: Diffusivity (D) obtained in simulations of heavy water for different functionals and in experiment, 
all at room temperature. Corresponding effective temperatures with respect to light water {T^^^^) and heavy 
water (r°f ). 

(xlO-5 cmVs) T^)f m T^^'ff (K) 



BLYP 0.20 240 (-20%) 250 (-17%) 

BLYP'' 0.15 237 (-21%) 245 (-18%) 

PBE"^ 0.16 238 (-21%) 246 (-18%) 

RPBE 1.50 283 ( -6%) 292 ( -3%) 

Exp (D2O)*' 1.87 

Exp (H2O)'' 2^50 




FIG. 2: Distribution of molecules with different coordinations. Room-temperature RPBE (circles, solid 
line) is compared with the TIP5P potential'*^ (diamonds, dashed line) at room temperature and with BLYP 
(squares, dotted line) at 345 K. 



as monitored by the heights of the first minimum and of the second maximum of gooi''')- The 
"moving- window" approach defined in Ref. [i^ has been used to compute the time averages. In the 
RPBE simulation windows of 5 ps were considered every 2 ps, whereas in the BLYP simulation the 
windows were of 7.5 ps every 2.5 ps, given the faster evolution of the former. Even if none of the 
trajectories has strictly evolved into a stationary situation, the figure reveals that the characteristic 
time of RPBE is shorter than the one of BLYP. However, from the tentative assimilation of the 
experimentally measured relaxation times^ to our equilibration times, one would expect for an 
effective temperature of 292 K a substantially shorted relaxation time. Surprisingly, it does not 
seem to be the case here, our equilibration time scale being rather on the 10 ps range or longer. 



IV. CONCLUSIONS 

We have performed DFT-based AIMD simulations of room-temperature liquid water, comparing 
the characteristics of the fluid obtained by BLYP and RPBE. The main conclusions are summarised 
as follows. 

(i) The RDFs for RPBE are closer to experiment than the corresponding ones for BLYP. RPBE 
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FIG. 3: Evolution of the difFusivity (a), the height of the first minimum of gooif) (b), and the concentration 
of coordination defects (c), in simulation 3 (RPBE) calculated in 5-ps- windows every 2 ps (solid line, circles) 
and in simulation 2 (BLYP) calculated in 7.5-ps-windows every 2.5 ps (dashed line, squares). 

water is less structured than BLYP, still slightly over-structured with respect to experiment. 

(ii) The diffusivity obtained for RPBE is still too low, but much closer to experiment than the 
value obtained for BLYP, with a factor of two underestimation instead of one order of magnitude. 

(iii) For room-temperature AIMD, the effective temperature (with respect to heavy water ex- 
perimental results) is ~292 K (-3%) for RPBE, while for BLYP it is ~250 K (-17%). In addition 
to the quantitative advantage of performing RPBE simulations there is a qualitative advantage in 
that neither of these RPBE simulations would be for supercooled water. 

(iv) The AIMD HB-network at AIMD temperature is extremely similar to the TIP5P network 
at the effective temperature. 

(v) The equilibration time-scale is around 10 ps or longer for room-temperature RPBE, as 
compared with more than 20 ps for BLYP. The fact that it is not substantially shorter is somehow 

• • 2(1 

surprismg. - 

Even if there are reasons to believe that the good results found for RPBE are not due to 
fundamental reasons, but rather to error cancellations, still, the possibility of having a well-behaved 
GGA water is appealing. We consider that it is worth further exploring the performance of RPBE 
in other wet systems, possibly starting from well characterised ion hydration shells. 
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